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I.  INTRODUCTION 


M-Layer  is  a  FORTRAN  program  for  computing  the  propagation  factor  of  an 
electromagnetic  (EM)  wave  in  a  stratified  atmosphere.  It  is  desirable  to  extend  the 
capability  of  this  program  to  include  a  layer  of  random  medium  representing  the  air- 
ocean  interface.  To  achieve  this  goal,  there  are  many  basic  theoretical  problems 
which  have  to  be  answered.  First  of  all,  the  effect  of  the  earth  curvature  in  this 
program  is  taken  care  of  through  the  classical  earth-fiattening  approximation  [Ref.  1], 
but  the  result  [Ref.  2]  does  not  agree  with  the  more  recent  diffraction  theory  of  Fock 
[Ref.  3]  near  the  surface  of  the  earth.  Then  there  is  the  question  about  the  better 
method  to  model  the  atmospheric  refractive  index  profile,  either  piecewise  linear  or 
quadratic,  to  be  resolved  by  a  new  earth-flattening  approximation  under  development 
at  NFS.  The  new  approximation  will  also  determine  the  functions  to  be  used  for  the 
representation  of  the  EM  fields  in  each  layer  through  uniform  asymptotic  theories. 
Within  some  proper  region,  these  new  functions  are  expected  to  reduce  to  the  Airy 
functions  utilized  by  M-Layer.  The  evolutionary  nature  of  this  effort  prompted  this 
review  to  improve  the  inner  workings  of  the  M-Layer  program.  In  particular,  the 
subroutines  to  search  for  the  modes  and  those  for  evaluating  the  Airy  functions  will 
remain  as  an  important  part  of  a  program  investigating  questions  about  EM  wave 
propagation  by  solving  the  related  boundary  value  problem. 
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It  can  never  be  overemphasized  that  a  boundary  value  problem  which  includes 
a  layer  of  random  medium  or  some  range  dependent  inhomogeneity,  set  up  according 
to  the  Maxwell  equations,  will  include  backscattering  in  its  solution.  This  is  in  sharp 
contrast  to  those  numerical  procedures  based  on  the  parabolic  approximation  to  the 
wave  equation  for  which  the  backscattering  is  completely  ignored. 

In  what  follows,  the  M-Layer  program  and  the  reasons  for  replacing  the 
extended  complex  numbers  with  their  complex  exponent  representations  are 
discussed,  together  with  some  other  problems  encountered  and  resolved  during  this 
investigation. 

A.  M-lAYER 

In  M-Layer,  the  index  of  refraction  of  the  atmosphere  is  assumed  to  be  height 
dependent  and  is  approximated  with  a  continuous  piecewise  linear  profile.  The 
classical  earth-flattening  approximation  is  utilized  to  allow  the  use  of  the  cylindrical 
coordinate  system  while  retaining  the  effect  of  the  curvature  of  the  earth.  This  is 
done  simply  by  substituting  the  index  of  refraaion  with  the  modified  index  of 
refraction,  which  also  has  a  piecewise  linear  profile  [Ref.  1]. 

The  source  of  the  EM  radiation  is  assumed  to  be  either  a  vertical  electric 
dipole  or  a  vertical  "magnetic  dipole",  with  the  latter  providing  an  approximation  to 
the  radiation  of  a  horizontal  electric  dipole.  The  dipole  is  located  along  the  positive 
z-axis  of  the  cylindrical  coordinate  system  while  the  origin  is  sitting  on  the  ground. 
The  x-y  plane  is  the  ‘flattened’  earth  surface.  After  carrying  out  the  Hankel 
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transfoim  along  the  radial  direction,  the  resulting  spectrum  of  the  Hertzian  dipole 
field  within  each  layer  of  a  linear  segment  of  the  modified  refractive  index  profile  is 
reduced  to  a  linear  combination  of  the  Airy  functions.  Specifically,  the  layers  are 
numbered  to  increase  with  height,  with  the  first  layer  being  the  one  above  the 
ground.  The  spectrum  of  the  Hertzian  dipole  field  is  proportional  to  the  product  of 
the  values,  at  the  transmitter  height  and  at  the  receiver  height  respectively,  of  the 
height-gain  function.  At  a  height  within  the  i-th  layer,  the  height-gain  function  is 
given  by  [Ref.  4]: 

where  p  is  the  radial  component  of  the  propagation  vector  and  is  also  the  spectral 
,  variable  of  the  Hankel  transform;  hence  it  is  the  same  throughout  all  layers.  It  is  a 
complex  variable  whose  imaginaiy  part  represents  the  radial  attenuation  rate  of  the 
spectral  component  of  the  Hertzian  dipole  field.  Under  the  classical  earth-flattening 
approximation,  the  spectrum  of  the  Hertzian  dipole  field  contains  a  discrete  portion 
and  a  branch  cut.  The  discrete  spectrum  gives  rise  to  the  creeping  wave  modes 
diffi’acted  by  the  earth  surface  and  the  dielectric  waveguide  modes  supported  by  the 
layered  atmosphere.  The  contribution  from  the  branch  cut  is  usually  negligible, 
especially  for  the  fit  in  the  shadow  of  the  earth.  The  M-Layer  program  locates  the 
discrete  spectrum  for  modes  having  a  radial  attenuation  rate  below  a  predetermined 
value.  Contributions  from  these  modes  determine  the  propagation  factor  of  the 
wave. 
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The  variable  qj  in  the  i-th  layer  is  a  dimensionless  linear  function  of  height  z 
with  the  free  space  wavenumber  k,  the  modified  index  of  refi’action  nij  at  the  lower 
boundary  z  =  z,-,  the  slope  of  the  modified  index  of  refraction  a,-/2  and  p  as 


parameters; 


2  0*^ 


(2) 


The  height  dependence  of  the  field  is  given  in  terms  of  the  functions  and 
^2(9,).  which  are  proportional  to  the  Airy  functions  and  Ai{-q^ 

respectively.  Of  these  two  functions,  at  a  height  so  large  that  9,-  is  large  and  positive, 
k-^fq^  represents  a  downward  going  wave  and  represents  an 

upward  going  wave.  The  coefficients^,-  andB,-  are  determined  by  the  conditions  on 
the  continuity  of  the  Hertzian  dipole  field  and  its  derivative  across  layer  boundaries 
and  by  the  normalization  condition  that  the  integral  of  the  square  of  the  height-gain 
function  over  all  height  equals  unity. 

To  fulfill  the  radiation  condition,  the  highest  layer  is  given  the  same  refractive 
index  as  the  free  space  above  it  and  only  the  outgoing  wave  is  allowed  within  this 
layer.  Below  the  ‘flattened’  earth  surface,  the  field  is  assumed  to  be  a  plane  wave 
propagating  downward.  Hence,  only  the  normalization  factors  are  required  in  the 
highest  layer  and  in  the  ground.  By  assigning  B,-  to  unity  in  the  highest  layer,  all  the 
coefficients  >4,-  andB,-  can  be  determined,  according  to  the  boundary  conditions,  to 
within  a  multiplicative  factor  forB,-.  This  multiplicative  factor  is  then  deduced  from 


the  normalization  condition.  This  procedure  can  also  be  carried  out  from  the  ground 
level  up.  That  these  coefficients  can  be  computed  either  from  the  highest  level  down 
or  from  the  lowest  level  up  is  a  result  of  the  fact  that  p  belongs  to  the  discrete 
spectrum  of  the  Hertzian  dipole  field.  Consequently,  agreement  between  these  two 
ways  of  evaluating  ihtA^  andS,-  coefficients  confirms  that  a  mode  has  been  located 
accurately. 

B.  EXTENDED  COMPLEX  NUMBER  REPRESENTATION 

The  discrete  spectrum  of  the  Hertzian  dipole  field  corresponds  to  the  zeroes 
of  the  modal  function  which  is  a  determinant  whose  elements  consist  of  and 
/c2(<7,)  at  the  layer  boundaries.  Numerically,  the  magnitude  of  this  modal  function 
causes  overflow  and  underflow  problems  as/l:j(<7,)  or  becomes  exponentially 
large  or  small  for  complex  qi  values.  In  the  M-Layer  program,  to  overcome  this 
problem,  a  complex  number  is  written  as  a  scaled  number,  which  is  complex, 
multiplied  by  a  scaling  factor  which  is  an  integer  power  of  e,  the  base  of  natural 
logarithm.  This  integer  is  chosen  so  that  the  greater  of  the  absolute  values  of  the 
real  part  and  the  imaginary  part  of  the  scaled  number  lies  within  e"^.  A  complex 
number  written  in  this  form  is  called  an  extended  complex  number.  Multiplication 
of  two  extended  complex  numbers  requires  summing  the  two  integer  exponents  in 
addition  to  carrying  out  the  regular  complex  multiplication  of  the  scaled  numbers. 
Addition  of  two  such  numbers  is  achieved  through  the  use  of  an  addition  subroutine: 
the  larger  scaling  factor  is  factored  out  of  both  addends  before  they  are  combined. 
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The  scaling  factor  is  adjusted  after  each  addition  and  after  a  sequence  of 
multiplications  to  make  sure  that  the  resulting  scaled  number  is  still  within  the 
desired  range.  Addition  is  troublesome  when  the  two  numbers  to  be  added  nearly 
cancel  each  other.  Under  this  circumstance,  the  scaling  factors  of  the  two  numbers 
are  identical  and  both  the  real  parts  and  the  imaginary  parts  of  the  scaled  numbers 
are  almost  equal  with  opposite  signs.  It  is  clear  that  the  real  part  and  the  imaginary 
part  of  the  sum  lose  their  accuracies  to  different  degrees;  hence  the  phase  angle  may 
incur  substantial  error.  To  remedy  this  situation,  interpolation  procedures  have  to 
be  devised. 

As  two  complex  numbers  come  close  to  cancel  each  other,  they  must  be  out  of 
phase  by  almost  180  degrees.  By  factoring  out  the  square  root  of  their  product 
instead  of  the  scale  factor,  the  resulting  addends  become  reciprocal  to  each  other, 
both  lying  within  an  identical  small  angle  to,  and  on  the  same  side  of,  the  imaginary 
axis.  They  are  close  to  the  unit  circle,  but  one  is  on  the  inside  and  the  other  is  on 
the  outside.  Taking  out  further  a  phase  factor  of  Tr/2  after  writing  the  addends  in 
their  exponential  forms,  the  exponents  become  small  numbers  for  which  a  Taylor 
series  expansion  of  the  exponential  function  converges  rapidly  and  can  be  used  for 
interpolating  the  sum  to  achieve  higher  accuracy.  Note  that  after  the  extra  phase 
factor  of  7r/2  is  removed  from  the  addends,  it  is  actually  the  difference  of  the 
resulting  two  reciprocals  which  is  computed.  This  procedure  effectively  picks  the 
direction  on  the  complex  plane  along  which  the  addends  are  almost  opposing  each 
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other  to  cany  out  their  cancellation.  The  resulting  sum  has  a  phase  angle  nearly 
perpendicular  to  this  chosen  direction. 

It  is  evident  that  the  representation  of  a  complex  number  by  its  complex 
exponent  of  base  e  provides  better  phase  accuracy  for  addition.  A  one-to-one 
correspondence  can  be  achieved  by  restricting  the  imaginary  part  of  this  exponent  to 
within  -JT  and  jt.  This  will  be  called  the  exponential  representation  or  the  complex 
exponent  representation  henceforth.  It  is  convenient  for  multiplication:  adding  the 
complex  exponents  of  the  two  factors  will  suffice.  Conversion  of  the  M-Layer 
program  from  the  extended  complex  number  to  the  complex  exponent  representation 
has  been  carried  out. 

C.  OTHER  REVISIONS 

As  better  precision  is  achieved,  problems  with  the  mode  search  procedure  and 
the  evaluation  of  the  and  5,-  coefficients  become  severe.  They  are  thoroughly 
investigated  and  resolved.  For  mode  search,  although  the  division  of  the  region  of 
interest  into  "contour  rectangles"  and  further  into  square  "meshes",  and  the  search 
pattern  to  move  around  the  sides  of  a  "contour  rectangle"  to  find  and  follow  "phase 
lines"  into  it  are  kept,  the  basic  assumption  of  Shellman  and  Morfitt  [Ref.  5]  that 
both  the  real  and  the  imaginary  parts  of  the  modal  function  are  linear  along  every 
edge  of  a  mesh  square  is  completely  abandoned.  For  the  evaluation  of  the  A,-  and 
Bj  coefficients,  the  "test  for  evanescence"  conditions  have  been  removed.  A  condition 
to  determine  whether  to  evaluate  the  coefficients  from  the  ground  level  up  or  from 
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the  top  level  down  has  been  fomulated  and  incorporated  into  the  program.  This 
accomplishment  leads  to  the  relaxation  of  mode  locating  accuracy  requirements 
which,  combined  with  the  improved  precision  of  the  revised  program,  makes  the  first 
order  Newton-Raphson  iteration  unnecessary.  The  specific  changes  in  the  program 
and  the  resulting  gains  in  speed,  accuracy  and  execution  stability  are  discussed  in  the 
following  chapters.  Suggestions  to  completely  revise  the  mode  search  protocol  to  do 
without  the  "contour  rectangles"  and  to  look  for  the  modes  according  to  their  range 
attenuation  rates  are  also  provided. 
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II.  PROGRAM  REVISIONS 


M-Layer  is  structured  into  three  parts:  setup,  mode  search  and  propagation 
factor  evaluation.  The  main  input  is  the  modified  refractive  index  values  at  specified 
heights  so  that  a  piecewise  linear  profile  can  be  constructed.  If  the  mode  locations 
for  the  particular  profile  are  available  fi’om  a  previous  run  of  the  program,  they  can 
also  be  included  in  the  input  and  the  mode  search  procedures  will  be  bypassed.  The 
various  ranges  and  transmitter  and  receiver  heights  for  which  propagation  factors  are 
desired  are  also  specified.  The  subroutine  WVGSTDIN  is  called  to  input  the 
information  from  an  ASCII  data  file.  The  program  then  computes  the  constants  to 
be  used  for  mode  search  and  propagation  factor  evaluation.  The  mode  search  is 
performed  with  the  subroutine  FNDMOD.  The  MODSUM  subroutine  is  then 
invoked  to  first  compute  the  A^-  and  coefficients  as  explained  in  the  Introduction, 
then  compute  the  propagation  factor  and  the  propagation  loss.  The  complete 
program  structure  is  given  in  Figures  1  and  2.  There  are  several  other  subroutines 
which  are  not  included  in  these  and  other  figures,  such  as  DHORIZ  for  computing 
the  horizon  distance  between  a  transmitter  and  a  receiver  for  reference  purpose; 
CHKMOD,  a  maintenance  routine  for  removing  zeroes  fi-om  reported  mode 
locations  by  older  versions  of  the  program;  or  A02H20,  a  routine  to  compute  the 
atmospheric  absorption  coefficient  due  to  oxygen  and  water  vapor.  They  will  not  be 
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Figure  1  Original  M-layer  subroutines  structure. 


Figure  2  Original  M<layer  subroutines  structure  (continued). 
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discussed  as  they  do  not  contribute  directly  to  the  main  purpose  of  this  program  of 
locating  the  modes  and  computing  the  propagation  factor. 

The  program  structure  has  been  altered  as  shown  in  Figures  3  and  4.  Since  the 
Aj  and  5,-  coefficients  have  to  be  evaluated  only  once,  they  are  now  obtained  through 
a  call  to  the  subroutine  ABCOEF  directly  from  the  main  program  right  after  the 
modes  are  located.  Several  subroutines  are  dropped  in  this  revision  for  various 
reasons:  The  subroutines  NORME  and  NORMRE  are  eliminated  because  they  are 
no  longer  needed  due  to  the  change  in  complex  number  representation;  the 
subroutines  NOMSHX,  FDFDTX  and  DXDETR  are  not  used  because  the  modes 
are  now  located  with  adequate  precision  without  further  iteration;  the  subroutine 
ADDX  is  not  listed  separately  because  it  is  called  only  once  and  has  been  reduced 
to  only  a  few  lines  which  are  placed  where  the  subroutine  is  called  in  the  original 
program.  On  the  other  hand,  changes  in  the  mode  search  algorithm  require  the 
addition  of  two  new  subroutines:  SURFO  is  a  modified  and  simpler  version  of  SURF; 
ROOTS  replaces  QUAD.  Due  to  the  change  in  complex  number  representation,  all 
subroutines  listed  below  FNDMOD  and  MODSUM  have  been  revised,  including 
their  input/output  lists.  But  except  for  SURFO  and  ROOTS,  the  utilities  of  these 
subroutines  are  the  same  as  those  of  the  original  ones.  Descriptions  of  these 
subroutines  can  be  found  in  the  report  by  Yeoh  [Ref.  4]. 

The  most  significant  changes  have  been  made  in  XCADD,  XCDATT  and 
XCDAIG  for  adopting  the  complex  exponent  representation  and  improving 
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Figure  4  New  M-layer  subroutines  structure  (continued). 
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computation  speed  and  accuracy;  in  FZEROX,  FINDFX,  ROOTS  and  SURFO  for 
stabilizing  and  simplifying  the  mode  search  algorithm;  and  in  ABCOEF  for 
implementing  the  criteria  to  determine  the  reliable  manner  for  evaluating  the  Af  and 
coefficients.  These  changes  are  discussed  in  the  sections  below.  The  source  code 
listings  of  the  completely  new  subroutines  XCADD  and  ROOTS  and  the  significantly 
revised  subroutines  FZEROX  and  ABCOEF,  which  are  compiled  with  Microsoft 
FORTRAN  version  5.00,  are  attached  as  Appendices  A  through  D.  Validation  of 
the  revised  program  has  been  carried  out  at  9.6  GHz  for  all  the  21  profiles  listed  in 
Yeoh  [Ref.  4]. 

A.  ADDITION  SUBROUTINE 

XCADD  is  the  subroutine  implementing  the  addition  of  complex  numbers 
under  the  representation  by  their  exponents.  Given  the  double  precision  complex 
numbers  Zj  and  Z2  as  the  exponents  of  the  addends,  this  subroutine  returns  the 
exponent  of  the  sum.  Since  a  double  precision  number  has  an  accuracy  of  53  bits, 
if  the  real  parts  of  Zj  and  Z2  differ  by  more  than  53  bits,  the  exponent  of  their  sum 
will  simply  be  the  one  of  the  greater  real  part.  When  cancellation  becomes  serious, 
the  square  root  of  the  addends  is  factored  out  first.  Then  the  four-term  Taylor  series 
expansions  of  the  resulting  reciprocals  are  summed.  Since  the  leading  term  of  the 
sum  of  the  Taylor  series  is  a  good  estimate  of  the  sum  of  the  reciprocals  and  the 
relative  error  of  the  four-term  Taylor  series  sum  is  proportional  to  the  fourth  order 
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of  this  leading  term,  the  threshold  for  invoking  this  interpolation  procedure  is  set  at 
the  highest  possible  value  of  2“^^  allowed  under  double  precision.  Experimenting 
with  this  procedure  shows  that  this  interpolation  improves  accuracy  as  long  as  the 
threshold  is  set  at  a  number  between  2”^  and  2~^^. 

B.  AIRY  FUNCTION  EVALUATION 

Similar  to  the  original  program,  the  evaluation  of  the  Airy  function  adopted  the 
algorithm  prescribed  by  Schulten,  et.  al.  [Ref.  6].  In  the  new  program,  changes  are 
made  to  follow  the  advice  of  Schulten,  et.  al.  concerning  the  region  within  which  a 
Taylor  series  expansion,  instead  of  the  faster  Gaussian  quadrature,  has  to  be  used  to 
achieve  double  precision  accuracy.  Other  changes  in  implementing  the  algorithm  are 
described  below. 

1.  XCDAIT 

Due  to  the  similarity  in  their  Taylor  series  coefficients,  the  Airy  function 
and  its  derivative  are  evaluated  within  a  single  loop.  The  relative  accuracy  of  the 
derivative  of  the  Airy  function  is  set  at  the  double  precision  limit  of  2“^. 

2.  XCDAIG 

Six  term  Gaussian  quadrature  is  used  for  evaluating  the  Airy  function  and 
its  derivative  outside  the  circle  of  radius  4.97  centered  at  (0.90, 2.80)  on  the  complex 
plane.  The  use  of  four-term  quadrature  outside  a  radius  of  15  from  the  origin 
suggested  by  Schulten,  et.  al.  is  not  adopted.  The  six-term  quadrature  in  this  range 
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retains  a  higher  accuracy  while  overall  speed  improvement  by  using  both  the  four- 
term  and  the  six-term  quadrature  appears  to  be  minimal. 

C.  MODE  LOCATING 

As  explained  in  the  Introduction,  the  modes  are  located  at  the  zeroes  of  the 
modal  function.  These  zeroes  are  located  on  the  upper  complex  plane.  Here  q^^ 
is  the  value  of  q^  on  the  earth’s  surface,  which,  according  to  Eq.(2)  of  Chapter  I,  is 
a  linear  function  of  p^.  For  a  horizontally  propagating  mode,  p/k  is  close  to  unity. 
The  maximum  range  attenuation  rate  specified  for  the  desired  modes,  which 
corresponds  to  a  limit  on  the  imaginary  part  of  p,  determines  approximately  the 
upper  bound  for  the  imaginary  part  of  the  q^  complex  plane  to  be  searched  for 
modes.  The  Shellman  and  Morffit  mode  search  procedure  first  divides  the  search 
region  horizontally  into  "contour  rectangles"  each  of  which  spans  160  meshes  along 
the  real  qii  direction.  A  mesh  is  a  square  whose  size  is  an  adjustable  parameter  of 
the  order  lO"'*  at  9.6  GHz  for  most  of  the  cases  considered  herein.  This  parameter 
is  determined  by  the  frequency  and  the  slope  of  the  modified  index  of  reflection  in 
the  lowest  layer  of  the  profile.  The  search  commences  at  the  top  left  comer  of  the 
"contour  rectangle"  whose  left  edge  has  a  real  coordinate  value  close  to  the 
difference  of  the  real  parts  of  the  q^j  values,  with  the  minimum  modified  index  of 
refraction  and  the  index  near  the  surface  substituted  into  Eq.(2)  of  Chapter  I.  After 
the  search  over  the  initial  rectangle  is  completed,  the  program  moves  to  search  the 
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next  rectangle  until  a  specified  maximum  number  of  modes  are  found  or  a  specified 
number  of  "contour  rectangles"  have  been  searched. 

The  search  for  zeroes  makes  use  of  the  fact  that  a  real  function  changes  sign 
when  it  crosses  a  simple  zero.  Since  a  zero  of  a  complex  valued  function  F(q)  is 
where  both  its  real  part  and  imaginaiy  part  vanish,  a  necessary  condition  for  a  point 
^  to  be  a  zero  is  that  it  is  on  the  intersection  of  two  curves  defined  by  Im{F(q)}  = 
0  and  Re{F(q)}  =  0.  The  program  searches  around  a  "contour  rectangle"  for  a  sign 
change  in  Im{F(q)}  across  an  edge  of  a  mesh  bordering  the  side  of  the  "contour 
rectangle"  to  determine  that  a  line  of  Im{F(q)}  =  0  has  been  encountered.  The 
search  then  follows  this  line  into  the  meshes  within  the  "contour  rectangle",  checking 
each  mesh  to  see  if  a  curve  Re{F(q)}  =  0  enters  the  mesh  under  investigation.  All 
these  steps  make  use  only  of  the  assumption  that  the  zeroes  of  the  modal  function 
are  simple.  Once  both  the  curve  Im{F(q)}  =  0  and  the  curve  Re{F(q)}  =  0  are 
determined  to  be  present  within  a  mesh,  the  location  of  their  possible  interception 
is  estimated.  An  algorithm  for  this  estimate  is  required. 

Shellman  and  Morffit  [Ref.  5]  introduced  a  further  assumption  that  the 
functions  Re{F(q)}  and  Im{F(q)}  are  both  linear  along  the  edges  of  a  mesh.  Based 
on  this  assumption,  they  try  to  estimate  the  locations  where  the  curve  Im{F(q)}  = 
0  enters  and  leaves  a  mesh  square,  and  the  location  of  if  a  curve  Re{F(q)}  =  0 
also  enters  the  same  mesh.  It  is  obvious  that  information  about  the  lo^tions  where 
the  curves  enter  and  leave  the  mesh  square  is  not  essential.  Furthermore,  in  the  18 
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m  duct  height  case,  the  scheme  causes  the  search  path  to  loop  around  four 
contiguous  meshes  until  the  search  is  broken  up  by  the  limit  on  the  number  of 
meshes  to  be  investigated.  Replacing  their  technique  requires  major  changes  in  the 
subroutines  involved.  A  new  subroutine  ROOTS  is  provided  to  estimate  the  location 
of  the  intersection  of  the  curves  Im{F(q)}  =  0  and  Re{F(q)}  =  0.  These  changes 
eliminate  the  looping  problem. 

Another  problem  is  encountered  in  the  40  m  duct  height  case  when  a  large 
number  of  zeroes  are  found  in  the  lower  half  complex  qj^  plane.  These  zeroes 
appear  to  belong  to  the  reflection  coefficient  on  the  wrong  sheet  of  the  branch  cut 
and  are  not  waveguide  modes.  This  happens  because  the  search  region  has  been 
extended  below  the  real  q^j  axis  to  avoid  the  singularity  in  SURF.  The  problem  with 
this  singularity  should  have  been  solved  within  SURF,  especially  because  it  occurs 
only  when  the  derivative  of  the  subroutine  output  variable  gamma  with  respect  to  q^ 
is  computed.  Since  this  derivative  is  not  needed  during  mode  search,  the  extension 
of  the  search  region  to  the  negative  q^  plane  is  unnecessary.  A  simplified  routine, 
SURFO,  is  introduced  which  is  exactly  the  same  as  SURF  except  that  it  does  not 
evaluate  the  derivative  of  gamma.  By  using  this  subroutine  instead  of  SURF,  the 
search  path  in  the  revised  program  does  not  avoid  the  real  and  the  imaginary  axes. 

1.  FNDMOD 

The  search  region  is  limited  to  the  upper  half  q^  plane.  All  the  modes 
found  are  ordered  according  to  their  range  attenuation  rates  before  those  numbered 
beyond  the  maximum  modes  allowed  are  abandoned. 
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2.  FZEROX 


Since  the  curve  Im{F(q)}  =  0  enters  into  a  mesh  square  through  an  edge, 
the  values  of  Im{F(q)}  must  change  sign  over  the  end  points  of  either  one  or  all 
three  other  edges.  When  there  is  only  one  other  edge  across  which  Im{F(q)} 
changes  sign  at  its  end  points,  it  is  the  edge  across  which  the  curve  Im{F(q)}  =  0 
exits  the  mesh  square.  Ambiguity  arises  when  all  edges  indicate  a  change  of  sign  at 
their  end  points.  When  this  occurs,  a  "right  turn  rule"  is  adopted  which  assumes  that 
the  curve  exits  the  edge  to  the  right  of  the  one  along  which  it  enters  the  mesh 
square.  Such  a  rule  avoids  the  retracing  of  the  search  path  when  the  mesh  square 
is  revisited  as  entering  this  same  mesh  square  from  the  left  side  of  an  edge  after 
exiting  from  its  right  side  requires  a  crossing  of  the  lm{F(q)}  =  0  curve,  which  is 
prohibited  under  the  simple  zero  assumption.  On  the  other  hand,  the  actual  curve 
may  have  turned  left  and  then  returns  to  this  mesh  square,  i.e.,  following  a  "left  turn 
rule.”  Under  such  a  scenario,  this  wrong  choice  would  have  left  a  segment  of  the 
curve  not  searched.  This  difficulty  has  not  been  observed  during  testing.  In  fact,  the 
ambiguous  situation  seldom  occurs.  Note  also  that,  as  remarked  above,  two  lines  of 
lm{F(q)}  =  0  do  not  cross  each  other  unless  a  higher  order  zero  is  present.  Hence, 
only  a  "right  turn  rule"  or  a  "left  turn  rule"  for  the  curve  to  exit  the  mesh  is  allowed. 
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Exiting  the  opposite  edge  demands  a  pair  of  crossing  Im{F(q)}  =  0  curves  within  the 
mesh  square.  This  violates  the  assumption  that  all  zeroes  are  simple.  Also  note  that, 
the  possibility  of  vanishing  Re{F(q)}  or  Im{F(q)}  values  at  the  comers  of  a  mesh 
square  is  eliminated  through  a  small  adjustment  in  FINDFX. 

3.  FINDFX 

Both  the  vertical  shift  away  from  the  real  axis  and  the  horizontal 
offset  away  from  the  imaginary  axis  are  uimecessaiy  and  have  been  removed  from 
this  routine.  Furthermore,  as  a  result  of  converting  to  the  complex  exponent 
representation,  the  sine  and  cosine  of  the  argument  of  the  modal  functions  are 
examined  for  sign  changes  in  FZEROX.  This  is  implemented  in  FINDFX  by 
including  the  cosine  and  sine  values  of  the  argument  of  the  modal  function  in  the 
output  list.  To  avoid  the  indeterminate  case  when  either  the  real  or  the  imaginary 
part  of  the  modal  function  becomes  zero  at  any  comer  of  a  mesh  square,  the 
argument  for  computing  the  cosine  and  sine  values  is  increased  by  2“^^  when  this 
occurs.  This  is  equivalent  to  a  consistent  small  distortion  of  the  particular  comer  of 
the  mesh  square.  This  will  not  cause  any  error  in  locating  the  zero  because  FINTDFX 
still  returns  separately  the  unmodified  exponent  of  the  value  of  the  modal  function. 

4.  ROOTS 

Assuming  that  the  modal  function  is  analytic  within  the  mesh,  this 
subroutine  utilizes  the  values  of  the  modal  function  at  the  four  comers  of  the  mesh 
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square  to  determine  the  Taylor  series  expansion  coefficients  of  the  modal  function 
to  the  third  order.  The  roots  of  this  cubic  polynomial  are  then  located  using 
Cardan’s  solution  by  radicals.  If  the  higher  order  coefficients  fall  below  machine 
resolution  for  a  root  within  the  mesh  square,  these  coefficients  are  regarded  as  zero 
and  the  order  of  the  polynomial  is  reduced  and  can  be  solved  more  expediently.  If 
the  function  is  determined  to  be  constant  over  the  mesh  square,  the  center  of  the 
square  is  taken  as  the  root  location. 

D.  EVALUATING  AND 

As  discussed  in  the  Introduction,  the  and  coefficients  can  be  evaluated 
either  from  the  top  level  down  or  from  the  lowest  level  up.  These  two  procedures 
are  simply  called  "integration  down"  and  "integration  up",  respectively,  in  the  original 
documentation  [Ref.  4].  The  location  of  a  mode  has  been  called  an  eigenvalue. 
That  the  results  of  "integration  down"  and  "integration  up"  agree  is  a  manifestation 
that  the  eigenvalue  is  located  accurately. 

The  subroutine  ABCOEF  evaluates  the  coefficients  A^  and  for  each  mode. 
If  the  range  attenuation  rate  for  a  mode  is  greater  than  0.1  dB/km,  the  coefficients 
are  evaluated  from  the  lowest  layer  up.  Otherwise,  it  is  evaluated  from  the  top  layer 
down.  It  is  obvious  that  such  a  rule  must  be  implemented  because  the  results  of 
"integration  up"  and  "integration  down"  do  not  agree  for  many  modes.  Efforts  are 
made  to  determine  the  cause  of  this  discrepancy  and  to  devise  a  means  to  resolve  it. 
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Investigation  reveals  that  inadequate  precision  in  the  location  of  the  modes  is 
one  source  of  the  problem.  Since  the  coefficients  depend  on  the  coefficients, 
while  the  Aj  coefficients  are  obtained  directly,  only  the  A^  coefficients  need  to  be 
examined.  The.<4j-  coefficients  of  the  six  modes  of  lowest  range  attenuation  rates  for 
all  21  profiles  except  the  one  without  evaporation  duct  are  computed  using 
eigenvalues  of  different  accuracy  controlled  by  the  first  order  Newton-Raphson 
iteration  method.  Table  1  shows  thtA^  coefficient  computed  with  the  new  program. 
They  are  arranged  from  the  top  layer  down.  In  the  i-th  layer,  the  Af  coefficient 
computed  by  "integration  downward"  depends  only  on>4,+;  in  the  layer  above  while 
that  computed  by  "integration  upward"  depends  only  on  A^^j  in  the  layer  below. 
Hence  in  each  layer,  the  coefficient  obtained  by  "integration  downward"  is  listed 
above  that  obtained  by  "integration  upward".  There  are  five  sets  of  A^  values  listed, 
with  the  magnitudes  given  in  powers  of  10,  and  the  phase  given  as  a  multiple  of  n. 
They  are  obtained  from  eigenvalues  of  decreasing  accuracy  -the  one  used  to  compute 
the  left  most  column  being  the  most  accurate.  The  first  set  is  computed  using  an 
eigenvalue  having  a  relative  accuracy  of  2"^.  The  second  set  uses  an  eigenvalue  with 
a  relative  accuracy  of  2“^.  The  relative  accuracy  of  the  eigenvalue  for  the  third  set 
is  2"^.  For  the  fourth  set,  the  first  order  Newton-Raphson  iteration  of  the  mode 
location  is  set  at  an  absolute  accuracy  of  0.03  of  the  mesh  size,  same  as  that  specified 
in  the  original  program.  The  eigenvalue  for  the  right  most  set  is  the  mode  location 
estimated  by  ROOTS  without  modification  by  the  Newton-Raphson  iteration.  It  is 
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TABLE  1.  IMPROVING  A,  ACCURACY  WITH  EIGENVALUE  18  M  DUCT 
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clear  that,  for  this  mode,  the  difference  between  these  two  methods  of  computing  the 
coefficients  becomes  negligible  as  the  accuracy  in  mode  location  increases.  For 
example,  in  the  8-th  layer,  the  magnitude  of  computed  by  integrating  downward 
changes  from  -1.9078  to  0.3482,  to  0.3460  to  03459,  which  agrees  with  the  result 
computed  by  integrating  upward.  The  phase  follows  the  same  trend  to  an  agreement 
within  0.001  n.  Table  2  shows  a  similar  set  of  output,  but  the  coefficients  fail  to  agree 
even  when  the  relative  accuracy  is  increased  to  2  Note  that  the  actual  difference 
in  both  the  real  part  and  the  imaginary  part  of  the  two  most  accurate  eigenvalues  is 
about  2  ■^.  Double  precision  accuracy  appears  to  be  insufficient  for  the  coefficients 
computed  with  these  two  methods  to  agree  for  all  modes.  Some  interesting  features 
can  be  observed  in  both  tables,  which  are  present  in  all  120  sets  of  values  computed. 
When  disagreement  is  present  in  one  set  of  Ai  coefficients,  such  as  those  in  either 
Table  1  or  Table  2,  the  change  toward  smaller  differences  with  improving  eigenvalue 
accuracy  occurs  mainly  in  one  way  of  computation,  but  not  both.  For  example,  in 
Table  1,  the  values  of  "integration  downward"  improve  with  better  eigenvalue 
accuracy,  while  those  computed  by  "integrating  upward"  change  little.  In  Table  2,  the 
results  of  "integration  downward"  are  the  ones  that  are  holding  steady  as  the  accuracy 
in  eigenvalue  improves.  Furthermore,  when  disagreement  occurs,  the  layer  in  which 
the  A^  coefficient  has  the  smallest  magnitude,  i.e.,  the  one  having  the  most  negative 
power  of  10,  divides  the  table  into  two  parts.  The  results  of  two  different  ways  of 
computation  agree  in  the  layers  above  this  one  if  they  disagree  in  those  below  it,  and 
vise  versa.  No  explanation  will  be  attempted.  Instead,  practical  rules  are  drawn  up 
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TABLE  2.  IMPROVING  A,  ACCURACY  WITH  EIGENVALUE  36  M  DUCT 
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to  take  advantage  of  these  facts.  In  Table  1,  the  process  of  "integration  upward"  goes 
through  the  troublesome  10-th  layer  and  produces  results  which  agree  with  the  results 
of  downward  integration  before  the  downward  process  goes  through  the  10-th  layer. 
On  the  other  hand,  the  downward  integration  is  tripped  up  going  across  the  10-th 
layer  and  produces  results  which  fail  to  agree  with  the  results  from  upward 
integration.  It  is  clear  that  the  results  from  upward  integration  are  the  correct  ones. 
This  conclusion  is  further  supported  by  the  fact  that  improving  the  accuracy  of  the 
eigenvalue  does  not  change  significantly  the  results  of  upward  integration.  Similar 
argument  leads  to  the  conclusion  that  in  Table  2,  the  results  of  downward  integration 
are  the  correct  values. 

It  can  be  concluded  from  the  above  observations  that  one  of  the  methods  of 
computing  the  coefficients  converges  to  the  correct  value  much  faster  than  the 
other.  It  is  also  found  that  this  method  of  faster  convergence  is  always  able  to  arrive 
at  the  correct  values  for  Af  for  all  the  cases  under  investigation. 

Table  3  lists  the  statistics  of  the  method  of  integration  which  yields  the  correct 
coefficients  for  each  of  the  120  modes  investigated.  The  differences  in  magnitudes 
and  phases  in  the  lowest  layer  and  in  the  layer  below  the  highest  are  also  listed. 
Since  for  most  of  the  cases,  when  disagreement  in  A^  values  occurs,  the  correct 
integration  is  upward  -this  is  used  as  the  default.  To  decide  that  downward 
integration  should  be  utilized,  the  following  steps  are  taken:  The  first  A^  value  of 
downward  integration  is  computed  and  compared  to  the  value  from  upward 
integration.  If  the  magnitudes  in  dB  disagree  by  less  than  0.02  dB,  their  phases  will 
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be  checked.  If  the  phases  differ  by  less  than  10~^)r,  the  agreement  is  deemed 
acceptable  and  the  Af  and  coefficients  computed  from  the  lowest  layer  up  are 
used.  Otherwise,  the  coefficients  are  re-evaluated  again  from  the  highest  layer  down. 

Once  the  correct  method  of  evaluating  the  and  coefficients  is  used,  the 
accuracy  of  the  mode  location  becomes  less  critical.  For  all  the  cases  investigated, 
the  A^  coefficients  obtained  from  mode  locations  estimated  with  or  without  the 
Newton-Raphson  first  order  iteration  differ  only  by  0.06  dB  in  magnitude  and 
0.0013  JT  in  phase  at  most.  In  fact,  few  cases  show  differences  more  than  0.002  dB 
and  0.0001  JT,  The  Newton-Raphson  iteration  is  not  needed.  Hence  the  subroutines 
NOMSHX,  FDFDTX  and  DXDETR  are  removed. 
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TABLE  3.  STATISTICS  FOR  EVALUATING  A,  COEFTICIENT 

A|A,|  m  Aarg(A,)/ii 

Duct  Mode  #  Evalnatiiig  Method 
hei^t 

Luju*  Layer 


up  I  down  bottmn  I  top-1  I  bottom  top-1 
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TABLE  3.  CONTINUED  1. 
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TABLE  3.  CONTINUED  2, 
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TABLE  3.  CONTINUED  4, 
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III.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  Performance 

This  revision  of  M-Layer  converts  the  extended  complex  number  representation 
of  an  exponentially  large  or  small  number  into  the  direct  representation  by  its 
complex  exponent.  The  accuracy  of  the  computation  has  been  improved  in  two  ways: 
First,  an  interpolation  algorithm  has  been  devised  when  severe  cancellation  of  the 
addends  is  detected.  Secondly,  accuracy  for  the  evaluation  of  the  Airy  function  has 
been  improved,  not  just  by  summing  the  Taylor  series  to  double  precision  resolution 
and  by  adopting  six-term  Gaussian  quadrature,  but  also  by  expanding  the  region 
within  which  the  more  expedient  Gaussian  quadrature  is  excluded  in  favor  of  the 
more  accurate,  but  time-consuming,  Taylor  series  summation.  The  improvement  in 
accuracy  is  most  easily  seen  from  Table  1. 

As  discussed  in  the  Introduction,  evaluating  the  Af  and  Bj  coefficients  either 
from  the  lowest  layer  up  (integration  up),  or  from  the  top  layer  down  (integration 
down),  must  result  in  the  same  values.  This  property  provides  a  consistency  check 
for  the  accuracy  of  the  computation.  For  the  six  modes  of  lowest  range  attenuation 
rates  of  the  20  profiles  of  different  duct  heights.  Table  1  lists  the  maximum 
difference  for  each  mode  which  shows  a  discrepancy  between  these  two  methods  of 
evaluating  the  Ai  coefficients.  For  each  profile,  the  maximum  value  in  magnitude 
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TABLE  1.  MAXIMUM  DIFFERENCE  IN  A,  COEFFICIENT  BETWEEN 
INTEGRATION  UP  AND  DOWN 


TABLE  1.  CONTINUED 


difference  in  dB  among  all  the  layers  is  listed  if  it  is  greater  than  2.  If  the  phases  of 
the  coefficients  deviate  more  than  O.Ijt  in  any  layer,  that  particular  mode  is  also 
singled  out.  The  location  of  the  mode  of  the  revised  program  is  within  a  relative 
accuracy  of  achieved  through  first  order  Newton-Raphson  iteration.  Even 
though  discrepancies  still  exist  when  the  duct  is  28  meters  or  higher,  it  is  clear  that 
the  revised  program  computes  more  accurately  than  the  original  one. 

For  the  cases  where  the  two  methods  of  evaluating  the  Af  and  B,-  coefficients 
disagree,  it  has  been  observed  that  one  of  the  methods  always  leads  to  values 
which  are  little  changed  when  the  accuracy  in  mode  location  is  varied,  while  the 
other  method  produces  A^  values  which  shift  toward  the  results  of  the  other  method 


as  the  accuracy  of  mode  location  improves.  Based  on  this  observation,  a  consistency 
check  is  implemented  into  the  program  to  identify  the  method  which  converges 
better.  For  the  120  cases  investigated,  when  this  method  of  faster  convergence  is 
used,  ihtA^  coefficients  obtained  from  mode  locations  estimated  with  or  without  the 
Newton-Raphson  first  order  iteration  differ  only  by  0.06  dB  in  magnitude  and 
0.0013  IT  in  phase  at  most.  In  fact,  few  cases  show  differences  more  than  0.002  dB 
and  0.0001  JT.  This  allowed  the  Newton-Raphson  iteration  to  be  removed  in  this 
revision. 

Table  2  compares  the  performance  between  the  original  and  the  revised 
programs.  The  time  spent  to  find  the  modes  has  been  reduced  by  an  average  of 
22.58%.  The  revised  program  can  always  produce  the  modes  found  by  the  original 
program.  Moreover,  the  mode  search  is  stable  for  the  new  program:  the  time  it 
requires  to  search  for  the  modes  is  about  the  same  for  similar  profiles.  The  sudden 
jumps  in  mode  search  time  for  the  24  m  and  the  40  m  cases,  which  indicate  troubles 
during  the  search,  no  longer  happen. 

With  the  proper  method  of  evaluating  the /I,  and  B,- coefficients  determined  by 
the  consistency  check,  the  output  of  the  revised  program  differs  fi’om  the  original 
program  in  some  cases.  The  most  serious  deviation  has  been  observed  for  the  38  m 
duct  height  case  as  shown  in  Tables  3  and  4.  For  example,  at  a  range  of  36.5  km 
with  the  transmitter  at  a  height  of  25  m  and  the  receiver  at  10  m,  the  coherent  path 
loss  is  175.93  dB  from  the  original  program,  and  is  167.90  dB  from  the  revised 
program. 
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TABLE  2  OVERALL  MODE  SEARCH  PERFORMANCE  COMPARISON 


DUCT 

HEIGHT 

ORIGINAL 

REVISED 

Time 

(meter) 

Time 

Modes 

Time 

Modes 

Improvement 

00 

0:00:37 

3 

0:00:35 

3 

5.40% 

02 

0:32:14 

9 

0:31:55 

9 

0.98% 

04 

1:14:12 

25 

1:05:04 

25 

12.31% 

06 

2:10:18 

53 

1:56:50 

53 

10.33% 

08 

0:35:58 

39 

0:29:25 

39 

18.21% 

10 

0:53:24 

59 

0:48:32 

61 

9.11% 

12 

1:09:40 

86 

1:01:44 

89 

11.39% 

14 

1:20:42 

94 

1:11:13 

97 

11.75% 

16 

1:54:35 

95 

1:18:07 

97 

31.82% 

18 

1:45:09 

100 

1:27:15 

104 

17.02% 

20 

1:46:19 

103 

1:34:20 

105 

11.27% 

22 

1:52:54 

105 

1:35:18 

106 

15.59% 

24 

3:42:59 

106 

1:46:47 

107 

52.11% 

26 

2:07:42 

106 

1:43:55 

108 

18.62% 

28 

2:00:05 

107 

1:44:59 

109 

12.57% 

30 

1:59:59 

107 

1:46:19 

108 

11.39% 

32 

1:55:29 

108 

1:42:58 

no 

10.84% 

34 

2:29:57 

109 

2:15:58 

111 

9.32% 

36 

2:31:40 

109 

2:17:20 

112 

9.45% 

38 

2:38:44 

110 

2:18:09 

111 

12.97% 

40 

5:41:17 

95 

2:39:39 

111 

53.22% 

Total 

40:23:54 

31:16:22 

22.58% 
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TABLE  3.  ORIGINAL  PROGRAM  38  M  DUCT  OUTPUT 

frequency  »  9600.0000  Rtiz 


renge 

Zt 

zr 

coherent 

incoherent 

coherent 

incoherent 

horizon 

(kn) 

(n) 

(M) 

node  um 

■ode  sue 

peth  loes 

peth  loss 

(km) 

(cB) 

(dB) 

(dB) 

(dB) 

27.3 

25.0 

4.0 

-15.30 

-15.62 

156.10 

156.43 

28.9 

27.3 

25.0 

6.0 

.62 

-2.35 

140.18 

143.16 

30.7 

27.3 

25.0 

8.0 

-1.11 

-4.21 

141.92 

145.01 

32.3 

27.3 

25.0 

10.0 

-27.26 

-12.66 

168.06 

153.46 

33.6 

36.5 

25.0 

4.0 

-16.94 

-16.62 

160.28 

159.96 

28.9 

36.5 

25.0 

6.0 

-.73 

-2.05 

1U.07 

145.39 

30.7 

36.5 

25.0 

8.0 

-2.21 

-3.72 

145.55 

147.06 

32.3 

36.5 

25.0 

10.0 

-32.59 

-14.29 

175.93 

157.64 

33.6 

45.8 

25.0 

4.0 

-19.89 

-16.96 

165.20 

162.26 

28.9 

4'!, 8 

25.0 

6.0 

-2.81 

-1.89 

148.11 

147.19 

30.7 

45.8 

25.0 

8.0 

-4.11 

-3.43 

149.41 

148.74 

32.3 

45.8 

25.0 

10.0 

-28.57 

-15.22 

173.88 

160.52 

33.6 

TABLE  4.  REVISED  PROGRAM  38  M  DUCT  OUTPUT 

frequency  *  9600.0000  mhz 


renge 

zt 

zr 

coherent 

incoherent 

coherent 

incoherent 

horizon 

(km) 

(m) 

(m) 

node  sum 

node  sun 

path  loss 

peth  loss 

(km) 

(dB) 

(dB) 

(dB) 

(dB) 

27.3 

25.0 

4.0 

-14.38 

-15.66 

155.18 

156.47 

28.9 

27.3 

25.0 

6.0 

.42 

-2.37 

140.39 

143.18 

30.7 

27.3 

25.0 

8.0 

-1.52 

-4.21 

142.33 

145.02 

32.3 

27.3 

25.0 

10.0 

-21.20 

-12.51 

162.01 

153.31 

33.6 

36.5 

25.0 

4.0 

-17.32 

-16.60 

160.66 

159.94 

28.9 

36.5 

25.0 

6.0 

-.48 

•2.08 

143.82 

145.42 

30.7 

36.5 

25.0 

8.0 

-1.62 

-3.73 

144.96 

147.07 

32.3 

36.5 

25.0 

10.0 

-24.56 

-14.04 

167.90 

157.38 

33.6 

45.8 

25.0 

4.0 

-20.26 

-16.93 

165.57 

162.23 

28.9 

45.8 

25.0 

6.0 

-3.14 

-1.93 

148.44 

147.23 

30.7 

45.8 

25.0 

8.0 

•4.62 

-3.46 

149.92 

148.76 

32.3 

45.8 

25.0 

10.0 

-25.40 

-14.90 

170.71 

160.21 

33.6 

37 


B.  Recommendations 


The  mode  search  protocol  of  this  program  needs  to  be  revised.  Since  the 
search  is  limited  by  the  number  of  modes  to  be  found  and  the  maximum  range 
attenuation  rate  accepted,  it  is  more  logical  to  begin  with  locating  the  mode  of  the 
lowest  attenuation,  and  then  proceed  to  look  for  the  next  one  in  the  order  of 
increasing  attenuation  rate.  Furthermore,  there  appears  to  be  only  a  single  ‘phase 
line’  of  vanishing  real  part  of  the  modal  function  on  which  all  the  modes  are  located. 
This  line  extends  from  lower  to  higher  range  attenuation  rates.  The  partition  of  the 
search  region  into  rectangles,  as  has  been  done  in  this  program,  tends  to  cut  the 
"phase  line"  into  segments  before  the  program  starts  to  search  for  the  end  points  of 
these  segments  and  then  follow  the  segments  in  different  directions.  It  is  clear  that 
a  better  way  for  mode  search  is  to  find  the  lower  end  of  the  single  "phase  line"  then 
follow  it  to  the  other  end. 
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APPENDIX  A;  SUBROUTINE  XCADD 


This  Appendix  lists  the  addition  subroutine  XCADD  which  returns  the  complex 
exponent  of  the  sum  when  the  complex  exponents  of  the  addends  are  given.  This  is 
a  complete  re-write  of  the  original  subroutine  of  the  same  name. 
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subroutfne  xcadd(zx,z1x,z2x) 


1 

2  c 

3  c  Given  zlx  and  z2x,  this  subroutine  adds  the  two  caif>lex  nuiters 

4  c  z1>exp(z1x)  and  z2>exp(z2x)  for  z>cxp(zx)  and  returns  zx. 

5  c 

6  c  inputs... 

7  c  zlxscooiplex  exponent  of  the  conptex  number  z1 

8  c  z2x>coinplex  exponent  of  the  coeptex  rxinber  z2 

9  c 

10  c  outputs... 

11  c  zx-conplex  exponent  of  the  cooplex  mmber  z 

12  c 

13  c  subroutines  called... 

14  c 

15  c********** 

16  inplicit  real*8  (a-h,o-z) 

17  co«plex*16  zx,z1x,z2x,zt1x,zt2x,clogzh,dsijn,czero,cerrx,cone,chpi 

18  paraineter(pi=3.141592653589793238462643d0,tijopi=2.d0*pi , 

19  hpis0.Sd0*pi,zeros0.d0,c16=1.d0/6.d0, 

20  *  bit14=1.d0/16384.d0,bit24sb{t14/1024.d0,ctot-bit14, 

21  >  4)is2259.cKI/4294967296.d0/4294967296.d0,hdpi>dpi/2.d0, 

22  ♦  e2n64=-3.742994775023704819d1.e2p27=-0.5d0*e2ii64, 

23  *  Chpis(0.d0,1.57079632679489661923132d0>.cone>(1.d0,0.d0), 

24  +  czero«<0.d0,0.d0),cerrx«<-3.742994775023704819d1,0.d0)> 

25  c  cerrx»e2ii64»-54*loa(2)aexponent  below  machine  accuracy 

26  dimension  ztnp(2),stiip(2) 

27  equivalence  (ztinp,closzh},(stiip,dsuni) 

28  c***** 

29  c  Replace  the  input  variables  with  a  local  variable  so  that 

30  c  equations  in  the  form  of  y=x*y  will  not  lead  to  confusion. 

31  c 

32  ztixszlx 

33  zt2x-z2x 

34  c 

35  Clogzh>0.5d0*(zt1x-zt2x) 

36  dxh«ztnp< 1 ) 

37  if(dxh  .It.  zero)  then 

38  zx«zt2x 

39  dxh*-dxh 

40  else 

41  zx«zt1x 

42  end  if 

43  c******** 

44  c  machine  accuracy  *  2**(*53) 

45  c  2**(27)«e**e2p27 
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46  c 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56  c 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67  c 

68  c 

69  c 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79  c 

80 


if  (dxh  .ge.  e2p27)  then 
return 
else 

zx>0 . 5d0*(  zt  1  X'i'Z  t2x  ) 
dsu»>cdexp<  c  I  ogzh ) 
dsuRBl  .d0/dsun»^un 
if  (cdabs(dsun)  .gt.  ctol)  then 
zxscdlog(dsum)>zx 
else 

Cencellation  is  serious.  latclogzhl  is  close  to  pi/2  or  *pi/2. 
yiBdnint(ztiRp(2)/tHopi  )*2.d0 
ztnp(2)>ztMp(2)*pi*yi 
dyi=d|pi*yi 

if  (ztMp(2)  .It.  zero)  then 
clogzh=*clogzh 
dyi=-dyi 
end  if 

ztinp(2}=(ztinp(2)-hpi  )-hdpi-dyi 

dsuRp2 . d0*c I ogzh* ( cone+c16*c I ogzh*c I ogzh ) 

if  (dsun  .eq.  czero)  then 

Note  that  a  complete  cancellation  of  two  nonzero  numbers  of 
order  one  is  considered  to  be  as  accurate  as  what  is  allowed 
by  the  machine  and  the  algorithm. 
zx>cerrx'^chpi>zx 
else 

dsumpcdlog(dsum) 

if  (stmpd)  .It.  e2m54}  stinp(1)=e2m54 
zx^dsun+chpi+zx 
end  if 
end  if 
return 
end  if 
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APPENDIX  B:  SUBROUTINE  FZEROX 


This  Appendix  includes  the  listing  of  the  subroutine  FZEROX  which  searches 
identifies  the  meshes  which  may  contain  modes  within  a  contour  rectangle.  The 
Shellman-Morffit  mode  locating  algorithm  has  been  completely  replaced. 


1  subroutine  fzerox(tleft, tright.tbot.ttop.tMshO, zeros, ni,nf) 

2  £**•** 

3  c  fzerox  is  a  routine  for  finding  the  zeroes  of  s  cosplex  function,  f, 

4  c  which  lie  within  a  specified  rectangular  region  of  the 

5  c  complex  q11  plww,  assuming  that  the  function  has  only 

6  c  simple  zeroes  over  this  rectangle. 

7  c 

8  c  parameters 

9  c  tleft  - 

10  c  tright- 

11  c  tbot  - 

12  c 

13  c  ttop  - 

14  c  tmesh  - 

15  c 

16  c 

17  c 

18  c  zeros  - 

19  c 

20  c  nf-ni  - 

21  c 

22  c  subroutines  calledd-* 

23  c  findfx 

24  c  roots 

25  c  nomshx 

26  e***** 

27  implicit  double  precision  (a*h,o-z) 

28  complex*16  f10,f01,f11,fxneu,fxold,fx00,fx10,fx01,fx11, 

29  czero, one, ci, sol, zeros 

30  parameter(czeros(0.d0,0.d0},one=(1.d0,0.d0),ci>(0.d0,1.d0)) 

31  $include:  'mlaparm. inc* 

*****  Begin  listing  of;  mlaparm. inc 

1  c 

2  c  include  file  to  define  the 

3  c  maximum  #  of  layers  (mxlayr) 

4  c  maxinuii  #  of  modes  (mxmode) 

5  c 

6  parameter  (mxlayr*35  ) 

7  parameter  (mxmode*127) 

End  listing  of;  mlaparm. inc 

32  dimension  kedge1(100),kedge2(100),keclge3(100},kedge4(100}, 

33  c  *  Ioc12r(mxmode),loc12i(mxmode),loc23r<mxmode),loc23i(mxinode), 

34  c  *  loc34r(mxfflode),loc34i(mxinode),loc41r(mxmode),loc41i(mxmode), 

35  *  sol(3),theta(2),zeros(2*mxmode*1) 

36  c 


specifying  the  search  rectangle; 
value  of  the  real  part  of  q11  at  the  left  edge, 
value  of  the  real  part  of  q11  at  the  right  edge, 
value  of  the  imaginary  part  of  q11  at  the  bottom  edge, 
(this  is  set  to  0.) 

value  of  the  imginary  part  of  q11  at  the  top  edge. 

set  equal  to  about  half  the  average  spacing  between 

zeroes  within  the  rectangle.  A  saialler  value  may  be  used 

as  a  safety  measure,  but  too  small  a  value  will  result 

in  excessively  long  run  time. 

output  list  of  (complex)  values  of  q11  at  which 

zeroes  are  found. 

the  number  of  zeroes  found 
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37  c 

38  comon  /tmccan/tnesh 

39  c***** 

40  c  Mxnsq  -  ■axinun  rxMber  of  Msh  aquares  allowed  on  any  one 

41  c  phase  line 

42  c  MBxnt  •  MxiMja  nuRter  of  tines  fzerox  will  reduce  tnesh 

43  c 

44  imxnscp3*MxO(int((ttop-tbot)/tiRshO).int((tright-tleft)/tiRshO)) 

45  maxnts? 

46  c***** 

47  tmesh  *  tinshO 

48  ntime  =  0 

49  go  to  7 

50  c 

51  5  tRieshstniesh/2.0d0 

52  ntime  =■  ntime«'1 

53  iflntime  .gt.  maxnt)  go  to  97 

54  c 

55  7  continue 

56 

57  c***** 

58  c  calculate  coordinates  of  rectangle  edges  in  tmesh  units 

59  c 

60  jit  s  idnint(tleft/tmesh'O.SdO) 

61  jrt  s  idnintftright/tmesh^^O.SdO) 

62  jtop  «  idnintCttop/tmesh^l.SdO) 

63  jbot  s  0 

64  c 

65  c  initialize  parameters  for  starting  search  at  tapper  left 

66  c  corner  of  search  rectangle 

67  c 

68  ki  s  Jtop 

69  kr  *  jit 

70  kedge  »  1 

71  call  f indfx(kr,ki,fxnew,xneu,yneu} 

72  nre1=0 

73  nre2s0 

74  nre330 

75  nre4=0 

76  knot12=0 

77  knot233:0 

78  knot34B0 

79  knot4l30 

80  nf=ni 

81  ni1*ni*1 
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90  to  15 


82 

83  c***** 

84  10  continue 

85  ifCnrzl  .It.  2)  go  to  15 

86  c  Hrite(16,2000)  nrzl 

87  go  to  5 

88  15  nrzl>0 

89  nrsqu  *  0 

90  20  fxoldsfxnew 

91  xol<^xneu 

92  yold^ynew 

93  go  to  (21,26,31.36). kedge 

94  £***•• 

95  c  search  along  left  edge  of  rectangle  for  changes  in  the 

96  c  sign  of  imagCf) 

97  c 

98  21  continue 

99  ifCki.eq. jbot)  then 

100  kedger2 

101  go  to  26 

102  end  if 

103  ki  *  ki-1 

104  call  f indfx(kr,ki ,fxnew,xneu,ynew) 

105  if  (yold*ynew  .gt.  O.dO)  go  to  20 

106  if(nrel.eq.O)  go  to  23 

107  c 

108  c  check  if  crossing  point  has  been  previously  found 

109  c 

110  do  22  k>=1,nre1 

111  if(ki.eq.kedge1(k))  go  to  20 

112  22  continue 

113  c 

114  c  follow  phase  line  through  rectangular  region 

115  c 

116  23  fx01>fxold 

117  fx01r>xold 

118  fx01i>yold 

1 19  f x00>f xnew 

120  fx00r>xnew 

121  fx00i«ynew 

122  li  ■  ki 

123  Ir  «  jit 

124  go  to  43 

125  c***** 

126  c  search  along  bottom  edge  of  rectangle  for  changes  in  the 
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127 

c 

sign  of  {iMg(f) 

128 

c 

129 

26 

continue 

130 

iftkr.eq.Jrt)  then 

131 

kedge«3 

132 

go  to  31 

133 

end  if 

134 

kr  ■  kr+1 

135 

call  findfxfkr.ki.fxnew.xneM.yneu) 

136 

if  (yold*yneM  .gt.  O.dO)  go  to  20 

137 

{f(nre2.eq.O)  go  to  28 

138 

c 

139 

c 

check  if  crossing  point  has  been  previously  found 

140 

c 

141 

do  27  k«1,nre2 

142 

if(kr.eq.kedge2(k))  go  to  20 

143 

27 

continue 

144 

c 

145 

c 

follow  phase  line  through  rectangular  region 

146 

c 

147 

28 

fxOOsfxold 

148 

fx00r«xold 

149 

fxOOi-yold 

ISO 

fx10»fxne« 

151 

fx10r«xnew 

152 

fx10i«ynew 

153 

li  •  jbot 

154 

Ir  «  kr-1 

155 

go  to  48 

156 

c**- 

157  c  search  along  right  edge  of  rectangle  for  sign  changes  in  iinsg(f>. 

158  c 

159  31  continue 

160  iffki.eq. jtop)  then 

161  kedge*4 

162  go  to  36 

163  end  if 

164  ki  «  ki«1 

165  call  findfx(kr,ki,fKnew,xneti,yneti) 

166  if  (yold*ynew  .gt.  O.dO)  go  to  20 

167  if(nre3.eq.O)  go  to  33 

168  c 

169  c  cheek  if  crossing  point  has  been  previously  fowxJ 

170  c 

171  do  32  k«1,nre3 
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172  if(k1.eq.kedge3(k})  go  to  20 

173  32  continue 

174  c 

175  c  follow  phase  line  through  rectangular  region 

176  c 

177  33  fx10«fxold 

178  fxIOr^xold 

179  fx10i-=yold 

180  fxiisfxnew 

181  fxllr^xnew 

182  fx11i>:ynew 

183  li  *  ki-1 

184  Ir  »  jrt-1 

185  go  to  53 

186  c***** 

187  c  search  along  top  edge  of  rectangle  for  sign  changes  in  iinag(f>. 

188  c 

189  36  continue 

190  iflkr.eq. jit)  go  to  80 

191  kr  »  kr-1 

192  call  findfx(kr,ki,fxnew,xnew,ynew) 

193  if  (yold*ynew  .gt.  O.dO)  go  to  20 

194  if(nre4.eq.O)  go  to  38 

195  c 

196  c  check  if  crossing  point  has  been  previously  found 

197  c 

198  do  37  k»1,nre4 

199  iflkr.eq. kedge4(k}}  go  to  20 

200  37  continue 

201  c 

202  c  follow  phase  line  through  rectangular  region 

203  c 

204  38  fxlixfxold 

205  fx11r*xold 

206  fx11i=yold 

207  fxOIsfxnew 

208  fx01r>xnew 

209  fx01i>ynew 

210  li  *  jtop-1 

211  Ir  ■  kr 

212  go  to  58 

213  c***** 

214  c  enter  stesh  square  from  left  side  or  exit  rectangle  at  right  edge. 

215 

216  41  lr»lr+1 
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217  If  (Ir  .1*.  Jrt'l)  go  to  42 

218  fwe3»nr^1 

219  kcdge3(nre3)«li'«'1 

220  go  to  10 

221  42  fx01>fx11 

222  fx01r»fx11r 

223  fx01<>fx11i 

224  fxOOsfxlO 

225  fx00r>fx10r 

226  fx00i«fx10i 

227  43  continue 

228  call  finelfx(lr+1,li+1,fx11,fx11r,fx11i) 

229  call  findfx<lr+1,li,fx10,fx10r,fx10i) 

230  c******* 

231  c  Determine  the  edge  of  exit  of  iffl(f)-0  from  current  mesh. 

232  edgeit=fx01i*fx11i 

233  edgeibsfx00i*fx10i 

234  if  (edgeib  .gt.  O.dO)  then 

235  c  Im(f)sO  goes  through  the  01  to  10  line. 

236  if  (edgeit  .gt.  O.dO)  then 

237  c  lm(f)=0  goes  through  the  10  to  11  edge  (edge  1). 

238  lout=1 

239  else 

240  c  lm(f)s0  goes  through  the  01  to  11  edge  (edge  2) 

241  (outs2 

242  end  if 

243  else 

244  c  lm(f)=0  goes  through  the  00  to  10  edge  (edge  4) 

245  lout=4 

246  if  (edgeit  .It.  O.dO)  then 

247  c  Im(f)30  also  runs  through  01  to  11  and  10  to  11  edges. 

248  c  Store  crossing  location  and  in/out  information. 

249  knot34>knot34*1 

250  c  loc34r(knot34)=lr 

251  c  loc34i(knot34)=li 

252  end  if 

253  end  if 

254  g******* 

255  go  to  60 

256  c***** 

257  c  enter  mesh  square  from  bottom  side  or  exit  rectangle  at  top  edge. 

258  46  li=li«1 

259  if  (li  .le.  jtop-1)  go  to  47 

260  nre4*nre4+1 

261  kedge4(nre4)*lr 


48 


262  go  to  10 

263  47  fx00«fx01 

264  fx00r>fx01r 

265  fx00i«fx01i 

266  fx10>fx11 

267  fx10r»fx11r 

268  fx10{=fx11i 

269  48  continue 

270  call  findfx{lr,li+1,fx01.fx01r,fx01i) 

271  call  findfx<lr+1,U+1.fx11,fx11r,fx11i) 

272  c******* 

273  c  Determine  the  edge  of  exit  of  im<f>*0  from  current  mesh. 

274  edgeil-fx00i*fx01i 

275  edgeir>fx10i*fx11i 

276  if  (edgeir  .gt.  O.dO)  then 

277  c  Im(f)s0  goes  through  the  00  to  11  line. 

278  if  (edgeil  .gt.  O.dO)  then 

279  c  lffl<f)=0  goes  through  the  01  to  11  edge  (edge  2) 

280  lout=2 

281  else 

282  c  Im(f)s0  goes  through  the  00  to  01  edge  (edge  3). 

283  loutx3 

284  end  if 

285  else 

286  c  lffl(f)s0  goes  through  the  10  to  11  edge  (edge  1) 

287  loutxl 

288  if  (edgeil  .It.  O.dO)  then 

289  c  Iin(f)x0  also  runs  through  00  to  01  and  01  to  11  edges. 

290  c  Store  crossing  location  and  in/out  information. 

291  knot41=knot41+1 

292  c  loc41r(knot41}xlr 

293  c  loc41i(knot41)=li 

294  end  if 

295  end  if 

296  c******* 

297  go  to  60 

298  c***** 

299  c  enter  mesh  square  from  right  side  or  exit  rectangle  at  left  edge. 

300 

301  51  lr=lr-1 

302  if  (Ir  .ge.  jit)  go  to  52 

303  nre1=nreU1 

304  kedge1(nre1)xli 

305  go  to  10 

306  52  fxll-fxOI 
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307  fx11r»fx01r 

308  fx11i*fx01i 

309  fxIO-fxOO 

310  fx10r»fx00r 

311  fx10i>fx00i 

312  S3  continue 

313  call  findfx(lr.li+1,fx01,fx01r,fx01i) 

314  call  findfxdr.li.fxOO.fxOOr.fxOOi) 

315  c******* 

316  c  Determine  the  edge  of  exit  of  im<f)»0  from  current  mesh, 

317  edgeit*fx01i*fx11i 

318  edgeib=fx00i*fx10i 

319  if  (edgeit  .gt.  O.dO)  then 

320  c  lin(f)E0  goes  through  the  01  to  10  line. 

321  if  (edgeib  .gt.  O.dO)  then 

322  c  lm(f)=0  goes  through  the  00  to  01  edge  (edge  3). 

323  lout=3 

324  else 

325  c  lm(f)*0  goes  through  the  00  to  10  edge  (edge  4) 

326  lout=4 

327  end  if 

328  else 

329  c  lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

330  lout>2 

331  if  (edgeib  .(t.  O.dO)  then 

332  c  Im(f)sO  also  runs  through  00  to  10  and  00  to  01  edges. 

333  c  Store  crossing  location  and  in/out  information. 

334  knot12=knot12+1 

335  c  loc12r(knot12)*lr 

336  c  loc12i(knot12)=li 

337  end  if 

338  end  if 

339  c******* 

340  go  to  60 

341  e***** 

342  c  enter  mesh  square  from  top  side  or  exit  rectangle  at  bottom  edge. 

343  56  li*li-1 

344  if  (li  .ge.  jbot)  go  to  57 

345  nre2=nre2'*’1 

346  kec^e2(nre2)>lr»1 

347  go  to  10 

348  57  fx01«fx00 

349  fxOIrsfxOOr 

350  fxOliofxOOi 

351  fxll-fxlO 
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352  fx11r«fx10r 

353  fx11i>fx10i 

354  58  continue 

355  call  findfx(lr,l{.fx00,fx00r,fx00i) 

356  call  findfx(lr»1.l{,fx10,fx10r,fx10{) 

357  c******* 

358  c  Determine  the  edge  of  exit  of  iBi<f>sO  froai  current  mesh. 

359  edgeil«fx00i*fx01i 

360  edgeir«fx10i*fx11i 

361  if  (edgeil  .gt.  O.dO)  then 

362  c  lM(f)>0  goes  through  the  00  to  11  line. 

363  if  (edgeir  .gt.  O.dO)  then 

364  c  Im<f)s0  goes  through  the  00  to  10  edge  (edge  4) 

365  lout>4 

366  else 

367  c  lffl(f)s0  goes  through  the  10  to  11  edge  (edge  1). 

368  lout=1 

369  end  if 

370  else 

371  c  lin(f)=0  goes  through  the  00  to  01  edge  (edge  3) 

372  lout=3 

373  if  (edgeir  .It.  O.dO)  then 

374  c  lffl(f)=0  also  runs  through  00  to  10  and  10  to  11  edges. 

375  c  Store  crossing  location  and  in/out  information. 

376  knot23»knot23+1 

377  c  loc23r(knot23)slr 

378  c  loc23i(knot23)-li 

379  end  if 

380  end  if 

381  c 

382  c******* 

383  uO  continue 

384  nrsqu=nrsqu+1 

385  if(nrsqu  .gt.  maxnsq)  go  to  95 

386  c****** 

387  c  Test  for  there  being  at  least  one  re(f)=0  line  entering  and 

388  c  leaving  the  aiesh  square. 

389  c 

390  if  ((fx00r*fx10r  .gt.  O.dO)  .and.  (fx01r*fx11r  .gt.  O.dO) 

391  +  .and.  (fx00r*fx01r  .gt.  O.dO))  go  to  (41,46,51,56)  lout 

392  c 

393  c  Computate  the  values  of  the  modal  function  at  the  corners  of  a 

394  c  a  mesh  square  to  determine  its  Taylor  series  to  the  3rd  order 

395  c  for  estimating  its  root  locations. 

396  c 
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397  c  fOO-one 

398  f10>cdexp(fx10-fx00ii-one 

399  f01>cdexp(fx01-fx00)-one 

400  f11«cdexp(fx11-fxOO)-one 

401  c 

403  c  write  (16,3001)  ni,nf,lr,l{,lcnot12,knot23,knot34,lcnot41 

404  c  3001  formt(/'  ni,  nf,  Ir,  li  and  knot12,  23,  34  and  43  before  ROOTS 

405  c  *•.'/,  2i6,2x,2i6,2x,4i6) 

406  c 

407  e***********  estiMte  locations  of  zeroes  by  radicals  ***************** 

408  c 

409  call  roots(f10,f01,f11,sol,nrsol) 

410  c 

411  do  63  n«1,nrsol 

412  ureal  «  dreaUsoKn)) 

413  uimag  >  dimagCsoUn)) 

414  if  (ureal  .It.  O.dO  .or.  ureal  .gt.  I.OdO)  go  to  63 

415  if  (uifliag  .It.  O.dO  .or.  uimag  .gt.  I.OdO)  go  to  63 

416  62  thetalDsdr^ureaD^tmesh 

417  theta(2)=(li‘»^uimag)*tmesh 

418  nf  *  nf*1 

419  zeros(nf  )>dciipl  x(  theta(  1 ) ,  theta(2) ) 

420  nrzl*nrzl+1 

421  63  continue 

422  c******************************************************************* 

423  c  write  (16,3002)  ni,nf,nrsol 

424  c  3002  formatf/'  out  of  ROOTS  at  63,  ni,  nf  and  #  of  roots  ',3i4) 

425  c******************************************************************* 

426  c  continue  following  the  phase  line 

427  go  to  (41,46,51,56)  lout 

428  c****** 

429  cc 

430  80  continue 

431  c 

432  return 

433  c***** 

434  95  continue 

435  wri tel 16,9500) 

436  wr i tel 16,4001 ) I r , I i , ni , nf , tmesh 

437  writel*  ,9500) 

438  4001  formatCgo  to  5  from  95  at  Ir,  li  *•, i6, i6, '  ni,  nf  =’,i6, 

439  ♦','fi6, ',  mesh  size  *',d14.6) 

440  go  to  5 

441  c***** 
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442  97  continue 

443  write(16,9700) 

444  write(16,4002)lr,li,ni,nf,tMesh 

445  write(*  ,9700) 

446  4002  forMt('90  to  5  from  97  at  Ir,  li  k',{6,',',{6,'  ni,  nf  >',i6, 

447  i6,',  Msh  size  *',d14.6,/'zeroe8  found  are  kept.') 

448  c  nf>ni 

449  c 

450  return 

451  c 

452  c****  forswt  stateawnts 

453  9500  foraiat(/5x, 'too  aiany  squares  on  same  phase  line  --  ', 

454  $  'reduce  tmesh  and  start  over') 

455  9700  fonnat(/5x, 'tawsh  has  been  reduced  but  problems  reamin  in', 

456  $  '  executing  fzerox') 

457  c 

458  end 
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APPENDIX  C:  SUBROUTINE  ROOTS 


This  Appendix  contains  the  listing  of  the  subroutine  ROOTS.  This  subroutine 
replaces  the  portion  of  the  subroutine  FZEROX  where  the  coefficients  of  a  quadratic 
equation  are  determined,  and  the  subroutine  QUAD  for  locating  the  zeroes  of  a 
quadratic  polynomial.  In  the  revised  subroutine  FZEROX,  the  roots  of  a  cubic 
polynomial  has  to  be  found.  This  subroutine  determines  these  zeroes  by  radicals. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 
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subroutine  roots  (f1,f2,f3,sol,nrsol) 

g******************************************************************** 

c  This  subroutine  finds  the  roots  of  a  third  order  polynomial  by 
c  radicals  when  the  values  of  this  polynomial  at  z«0,  z«1,  z«i  and 
c  z«1«i  are  given  as  f0=1,  fUfO,  f2*f0  and  f3«^f0  respectively, 
c  Note  that  this  algorithm  takes  cubic  roots  of  two  complex  nuabers 
c  (hence  the  name  'solution  by  radicals')  and  use  their  linear 
c  combinations  as  the  roots  of  a  third  order  polynomial. 

implicit  real*8  (a-h,  o-z) 

cofflplex*16  f 1, f2, f3, zero, one, ci,ep14,effl14,ep23,em23, 

fa,fb,fc,fd,fa1,fa2,fa3,fa1s,p,q,delt,z,zffl,u,v,sol 
parameter  (xbit52=52.d0*0.69314718055994531d0,thrd=1.d0/3.d0, 
bit50=1 .d0/33554432.d0/33554432.d0,bit51=bit50/2.d0, 
♦  bit52=bit51/2.d0,tol=0.001d0, 

zero=(0.d0,0.d0),one=(1.d0,0.d0),cis(0.d0,1.d0), 

+  ep14=(0.5d0,0.5d0),effl14=(0.5d0,-0.5d0), 

ep23=(-0.5d0,0.86602540378443864675d0), 
em23=(-0.5d0,-0.86602540378443864675d0)) 
dimension  sol(*) 
fa=one 

fb=(f2-ci*f1+em14*f3) 
fc=<(ep14*-one)*f1-(em14+one)*f2+ci*f3) 
fd=(em14*<f2-f1)-ep14*f3) 
if  (cdabs(fb)  .le.  bit50)  fb^zero 
if  (cdabs(fc)  .le.  bit51)  fc^zero 
if  (cdabs(fd)  .le.  bit52)  fd^zero 
if  (fd  .ne.  zero)  then 
fa1=(-thrd)*fc/fd 
fa2=fb/fd 
fa3=fa/fd 
fa1s=fa1*f8l 
p=thrd*fa2-fa1s 

q=0.5d0*<fa3+fa1*fa2)-fa1*fa1s 
if  (p  .eq.  zero)  then 
if  (q.  eq.  zero)  then 
nrsol=1 
sol(1)=fa1 
return 
else 

nr8ol=3 

u=((-2.d0)*q)**thrd 
sol(1)*u*f8l 
sol(2)=ep23*u»fa1 
sol(3)»em23*u»fa1 
return 
end  if 
else 
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49 

if  (q.  aq.  zero)  then 

50 

nrcol>3 

51 

aoU1)>fa1 

52 

u>cdsqrt( ( ■3.d0)*p) 

53 

eol(2)>fe1'HJ 

54 

soU3)*fa1-u 

55 

return 

54 

else 

57 

v»p/q 

58 

z«p*v*v 

59 

ebez>>cdabe(2) 

60 

if  (ebsz  .It.  tol)  then 

61 

zms-z 

62 

f n>di nt ( 1 .dO-xbi t52/dlog(absz) ) 

63 

lsstn«idint(fn)-1 

64 

drm^fn-O.SdO 

65 

dnd=fn+1.0d0 

66 

delt=one 

67 

do  100  nt«1,lastn 

68 

dnnsdnn-l.dO 

69 

dnd=dnd-1.d0 

70 

del  t*<dnn/dnd)*del  t*zint-one 

71  100 

continue 

72 

delt=<0.5d0*delt/q)**thrd 

73 

u=p*delt 

74 

V3-1.d0/delt 

75 

else 

76 

de 1 t «cdsqr t < one+ z ) - one 

77 

u=<q*delt)**thrd 

78 

ys.p/u 

79 

end  if 

80 

nrsolO 

81 

sol<1)=u+v+fa1 

82 

sol (2)=ep23*u*em23*v+f a1 

83 

sol  (3)®eni23*u*-ep23*v*f  a1 

84 

return 

85 

end  if 

86 

end  if 

87 

else  if  (fc  .ne.  zero)  then 

88 

if  <fb  .eq.  zero)  then 

89 

if  (fa  .eq.  zero)  then 

90 

nrsol*1 

91 

sol(1)>zero 

92 

return 

93 

else 

94 

nrsol>2 

95 

z»cdsqrt(-f8/fc) 

96 

sol(1)>z 
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97 

sol(2)=-z 

98 

return 

99 

end  if 

100 

else 

101 

fa1«0.5d0*fb/fc 

102 

fa2-fa/fc 

103 

ZBfa2/fa1/fa1 

104 

afa6Z3cdafas(z) 

105 

if  (absz  .It.  tol)  then 

106 

fn>dint(1  .dO-xbi  t52/dl09(ab6Z)) 

107 

lastnBidint(fn)-1 

108 

dnn>fn-0.5d0 

109 

dnd>fn»1.0d0 

110 

deltsone 

111 

do  200  ntsl.lastn 

112 

dmadnn-l.dO 

113 

chd=dnd-1.d0 

114 

del  ts(dnn/dnd)*del  t*Z'*'one 

115  200 

continue 

116 

delt=-0.5d0*delt/fa1 

117 

nrsol-2 

118 

sol(1)sfa2*delt 

119 

sol(2)-1.d0/delt 

120 

return 

121 

else 

122 

de 1 tscdsqrt (one- z ) 

123 

nrsols2 

124 

sol<1)=-fa1*(one-delt) 

125 

sol<2)=-fa1*{one+delt) 

126 

return 

127 

end  if 

128 

end  if 

129 

else  if  (fb  .ne.  zero)  then 

130 

nrsol=1 

131 

sol(1)=-fa/fb 

132 

return 

133 

else 

134 

nrsol=1 

135 

sol(1)zep14 

136 

return 

137 

end  if 

138 

end 
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APPENDIX  D:  SUBROIHINE  ABCOEF 


This  Appendix  contains  the  listing  of  the  subroutine  ABCOEF.  The  consistency 
self-checking  procedure  has  been  implemented  to  determine  the  correct  method  to 
evaluate  the  and  B,-  coefficients. 
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subroutine  abcoef(zero,in) 


1 

2  c’ 

3  c  For  each  node  m,  this  suboutine  calculates  A-B  coefficients  in 

4  c  all  layers  for  combining  two  linearly  independent  solutions  of 

5  c  Stokes'  equation  to  fona  the  height  gain  function: 

6  c 

7  c  height  gain>exp(bcoefx(l,«i))*(k1*exp(acoefx(l,a))«k2) 

8  c 

9  c  where  k1  and  k2  are  tuo  independent  solutions  to  Stokes' 

10  c  equation.  In  the  top  layer  (i.e.  nzlayr)  the  height  gain  is: 

11  c 

12  c  height  gainsexp(bcoefx(l,M))*h2 

13 

14  c  where  h2  is  a  solution  to  the  Stokes'  equation  associated 

15  c  with  outgoing  energy  flow.  Here  k1  and  k2  are  proportional 

16  c  to  the  k1  and  k2  used  by  Marcus  and  the  h2  is  proportional 

17  c  to  a  modified  Hankie  function  of  order  1/3. 

18 

19  c  inputs... 

20  c  zero-an  eigenvalue  in  q11  space 

21 

22  c  outputs... 

23  c  acoefx-two  dimensional  array  of  complex  exponents 

24  c  coefficients  used  to  combine  two  linearly 

25  c  independent  solutions  of  stokes'  equation 

26  c  bcoefx'two  dimensional  array  of  complex  exponents 

27  c  coefficients  used  for  normalizing  the  height  gains 

28 

29  c  note:  ecoefx  and  bcoefx  are  passed  by  the 

30  c  common  block  /pap2/ 

31 

32  c  subroutines  called... 

33  c  xcdai 

34  c  xcadd 

35 

36  c  common  block  areas... 

37  c  comi 

38  c  coffl2 

39  c  papi 

40  c  pap2 

41  c***** 

42 

43  implicit  real*8(a-h,o-z) 

44  cofflplex*16  acoefx, bcoefx, Cqij,h2xq1,(ti2xq1,h2xq2,dh2xq2.k1xq1, 

45  $  dk1xq1,k1xq2,dk1xq2,k2xq1,dk2xq1,k2xq2,dk2xq2,h2dk1x, 
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46  t  clt2lc1x,h2dk2x,dh2l(2x,nuMX,(lenax.niab)(,denbK,int1x,int2x, 

47  S  hyx,dh>-x,k1dhyx,dk1hyx,dk2hyx.k2dhyx.9«M,dBai«dq,<, 

48  *  koa123,rtsuiix,zero,q1,q2,«iiix,surfno,dqij,dq{jdz,sqn9, 

49  S  dnuit)x,dhux,dhlx,e13x,cneg,cldqzl,cldq»,c<gaM,koaMv,tthd, 

50  *  tacoef.dacoef 

51 

52  para«eter(doiini«1.d-3,do*<nr»1.d-3/0.4342944819032518d0, 

53  «  pi33.141592653589793238462643d0, 

54  ♦  i«(0.0d0,1.0d0>,tthd«(2.d0/3.d0)*i, 

55  cne9>(0.0d0,3.141592653589793238462643d0),e13xxcneg/3.d0) 

56  c***** 

57  c  mxlayr^iMxiiiui  rxnber  of  layers  alloweo 

58  c  inxmodesfflaxiiixin  number  of  inodes  allowed 

59 

60  c 

61  c  use  include  file  for  parameters  of 

62  c  use  include  file  for  parameters  of 

63  c  mxlayr  max  #  layers 

64  c  mxmode  max  #  modes 

65  c 

66  Sinclude:  'mlaparm.inc' 

*****  Begin  listing  of:  mlaparm.inc 

1  c 

2  c  include  file  to  define  the 

3  c  maxinun  #  of  layers  (mxlayr) 

4  c  maximum  §  of  modes  (mxmode) 

5  c 

6  parameter  (mxlayr=3S  ) 

7  parameter  (mxmode=127) 

*****  End  listing  of:  mlaparm.inc 

67  c 

68  c 

69  c***** 

70  c  acoefx-two  dimensional  complex  array  used  for  combining  two 

71  c  independent  solutions  to  stokes*  equation 

72  c  bcoefx-two  dimensional  complex  array  used  for  normalizing  height 

73  c  gain 

74  c  cqij-two  dimensional  array  containing  coefficients  for  evaluating 

75  c  qij  in  terms  of  q11 

76  c  dqij-array  containing  coefficients  for  evaluating  qij  in  terms  of 

77  c  q11 

78  c  dqijdz-array  containing  derivatives  of  qi(z)  in  the  different 

79  c  layers 

80  c  zi-array  containing  input  hesights  for  the  modified  refractivity 

81 
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82  dimension  acoefxdwlayr.Mxmode), 

83  %  bcoefx(mxlayr,mxmode), 

84  S  dqi  j(iiixlayr),cqi  j(mxlayr.2).dqi  jdzCmxlayD.zKmxlayr+l) 

85  c***** 

86 

87  caamon  /coai1/freq,Haveno,sqng 

88  coninon  /com2/cq1  j,dqi j,dqi jdz.nzlayr 

89  coomon  /pap1/nrmode,koa123,8urfno,zi 

90  common  /pap2/acoefx,bcoefx 

91 

92  c***** 

93  c  check  for  single  layer 

94  c 

95  c  set  e  complex  variable  koawav«*i*koa123/(waveno*Haveno)  to 

96  c  avoid  repeating  computations 

97 

98  koaMav=-i*koa123/(uaveno*waveno) 

99 

100  if(nzlayr  .eq.  Dthen 

101  q1=cqi j(1,1)*zero*dqi j<1) 

102  call  surf(q1,gannia,dgaffldq) 

103  call  Xcdai('q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 

104  dh2xq1>dh2xqUe13x 

105  inti  x»cd  I  og(  koawav*dgaffldq-q1  /'iqi  jdz(  1 )  >*2 .  OdO*h2xq1 

106  int2x>2.0d0*dh2xq1 -cdlog( -dqi Jdz(1 ) ) 

107  call  xcadd(sunx,int1x,int2x) 

108  rtsumx=0.5d0*sijnx 

109  bcoefx(1,m)s*rtsiinx 

110  return 

111  end  if 

112 

113  cldqzl=cdlog(-dqi jdz(l)) 

114 

115  c  if  I  equals  one  then  initialize  cunulants  and  caculate  a's  and 

116  c  b's  in  bottom  layer  using  ground  boundary  conditions. 

117 

118  q1=cqi j(1,1)+zero*dqi j<1) 

119  call  Xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 

120  dk2xq1>dk2xq1'»cneg 

121  dk1xq1^dk1xq1-e13x 

122  call  surf (q1, gamma, dgamdq) 

123  cigama>cdlog(i*gamma) 

124  call  xcadd(numax.cldqzl*cneg+dk2xq1.cigame^cneg>k2xq1} 

125  call  xcadd(denax,cigaffla«k1xq1,cldqzUdk1xq1) 

126  acoefx(1,m)>numax-denax 
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127 

128 

129 

130  c 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166  c 

167 

168 

169 

170 

171 


call  xcadd(clenbx,k2xq1,acoefx(1,i*)-i'lc1xq1) 
bcoef  x(  1  ,in)>*denbx 

calculate  contributions  to  nonaalizing  integrals. 

ca  1 1  xcadd(  hyx,  k2xq1 ,  acoef  x(  1  .aD-^klxql ) 
hyx=bcoef  x(  1  .ml-t-hyx 

ca  1 1  xcadd(dhyx ,  dk2xq1 ,  acoef  x(  1  ,ffl)‘»dk1  xql ) 
dhyx-bcoef  x(  1  .■D'fdhyx 

i  nt  1  x-cdl  og(  koauav*dgaflidq*q1  /dqi  jdz(  1 )  )«^2 . 0d0*hyx 
i nt2xs2 . 0d0*dhyx- c Idqz I 
call  xcadd(siinx,int1x,int2x) 

do  9010  l»2,nzlayr-1 
lin1=l-1 

cldqzl=cdlog(-dqi jdz(l)) 
c  I  dqztn=cd  I  og  (  dq  i  j  dz  ( I  ml )  ) 
q1=cqi j  < 1 , 1 )+zero*dqi j ( I ) 

call  xcdai (-q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1 =dk2xq1 +cneg 

dk1xq1=dk1xq1-e13x 

q2=cqi j ( Iml ,2)+zero*dqi j ( Iml ) 

call  Xcdai(-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2'»cneg 

dk1xq2:cdk1xq2-e13x 

call  xcadd< hyx , k2xq2 , acoef x( Iml , m)'<'k1  xq2 ) 

call  xcadd(dhyx,dk2xq2,acoefx( Iml ,m}+dk1xq2) 

kldhyx^klxqUdhyx 

dk  1 hyx=dk 1 xql+hyx 

dk2hyx=dk2xqUhyx 

k2dhyx=k2xq1*dhyx 

cal  I  xcaddldenax.cldqzKH'kldhyx.cldqzUdklhyx) 

cal  I  xcadd(numax,cldqzl-cneg+dk2hyx,ctdqzm4’cneg+k2dhyx} 

acoef x( I ,m}=numax-denax 

cal  I  xcadd(denbx,k2xq1  ,acoefx(  I  .nD-^klxql } 

rMii>x=bcoefx(  Iml  ,m)'*'hyx 

dnumbxsbcoefxl Iml ,m)+dhyx 

bcoefxl  I  ,m)=nutnbx*denbx 

calculate  contribution  to  normalizing  integrals. 

intlx^cdlogC  -ql/dqi  jdz(  I  )'*’q2/dqi  jdz(  Im1  ))'»2.0d0*nLirbx 
call  xcadd(sunx,suiix,int1x) 
call  xcadd(dhux,dk2xq1  .acoefxd.ml^dklxql ) 
dhux=bcoef  x(  I  ,m}'«dhux 
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172 

i  nt  1  x>2 . 0d0*dnu«fax  -  c  1  dqzin 

173 

i nt2xs2. 0d0*dhux-c Idqz 1 

174 

call  xcadd<sumx,8iiiix,int1x) 

175 

call  xcadd(suiix,sijiix,int2x) 

176 

9010  continue 

177 

% 

178 

179 

c  if  1  equals  nzlayer,  calculate  a's  and  b's  using  outgoing 

180 

c  wave  in  top  layer. 

181 

182 

nzni1»nzlayr-1 

183 

q1=cqi j(nzlayr,1)*zero*dqi j(nzlayr> 

184 

cal  1  xcdai ( -q1 , k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

185 

dh2xq1  sdh2xq1  ■•’e13x 

186 

q2scqi  j  (nzmi  ,2)'*'zero*dqi  jlnzml ) 

187 

call  xcdai (-q2,k2xq2,dk2xq2,k1xq2,dk1xq2.h2xq2,dh2xq2) 

188 

dk2xq2=dk2xq2-i'Cneg 

189 

dk1xq2sdk1xq2-e13x 

190 

cal  1  xcadd(hyx.k2xq2,acoefx(nzffl1.in)-»k1xq2) 

191 

mjnbx^bcoef  x(  nz  1  ayr- 1 ,  ml^hyx 

192 

bcoef  x(nz  Iayr,ni)=nijik>x-h2xq1 

193 

194 

c  calculate  contribution  to  cunulants. 

195 

196 

i nt 1 x=cdlog< -ql/dqi jdzlnzlayr )+q2/dqi jdzlnzml ) )♦ 

197 

$  2.0d0*numbx 

198 

call  xcadd(siinx,sijnx, intlx) 

199 

call  xcadd( dhyx , dk2xq2 , acoef x( nzml , in)+dk  1  xq2 ) 

200 

dnumbx=bcoefx(nzin1  ,ni)+<#iyx 

201 

int1x=2.0d0*dnuit)x-cdlog(dqi  jdz(nzm1  >) 

202 

call  xcadd(sumx,sunx, intlx) 

203 

dhuxsbcoef  X  ( nz  1  ayr ,  in)+dh2xq1 

204 

int2x=2.0d0*dhux-cdlog(-dqi jdzfnzlayr}) 

205 

call  xcadd(sunix,8utnx,int2x) 

206 

207 

c  renormalize  b's  so  that  height  gain  integral  equals  unity. 

208 

209 

rtsumxs . 5d0*sunx 

210 

do  9000  ll^l.nzlayr 

211 

bcoefxl 1 1 ,m)sbcoefx( 1 1 ,m)-rtsunx 

212 

9000  continue 

213 

214 

215 

% 

216 
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217  l*nzlayr 

218  l«1»l-1 

219  cldqznBC<ilo9(dqi  jdz(lml)) 

220  clclqzl>cdlog(*dqi  jdz(l)) 

221 

222  c  calculate  q  and  associated  quantities  at  bottoai  of  layer  I 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233  c****» 

234  c 

235  c 

236  c***** 

237 

238 

239 

240 

241 

242 

243 

244 

245  c  If  In  the  nzlayr-1  layer  the  magnitudes  of  A  coefficients  from 

246  c  Integration  up  and  down  differ  by  less  than  0.02  dB  and  their 

247  c  phases  differ  by  less  than  O.OOlpi,  the  A  and  B  coefficients 

248  c  obtained  from  integration  up  will  be  accepted. 

249 

250  tacoef=rMnaX'denax 

251  dacoef=tacoef 'acoefx(lffl1,m) 

252  dlfr^bs(dreal(dacoef )) 

253  if  (difr  .It.  downr)  then 

254  dif isdlmagldacoef )/pi 

255  difi=dabs(difi-dnint(difi/2.d0)»2.d0) 

256  if  (difl  .It.  douni)  return 

257  end  if 

258 

259  acoefx(lm1,m}>t.scoef 

260  call  xcadd(denbx,lc2xq2,acoefx(lm1,m)«k1xq2} 

261  bcoefx(lm1,m)3h2xq1-denbx 


qlscql j< I , 1 )*zero*dql j< I ) 

call  xcdai ( -q1 , k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,(ti2xq1 ) 
dh2xq1 3dh2xq1 '•'e13x 

q2=cql  j  ( Iml  ,2)«zero*dql  j  ( lail ) 

call  Xcdai(-q2,k2xq2,dk2xq2,k1xq2.dk1xq2,h2xq2,dh2xq2) 

dk2xq2sdk2xq2'«'Cneg 

dk1xq2=dk1xq2-e13x 

Caculate  acoefx( Iml ,m),bcoefx( tml ,m) 
and  cumulants  using  outgoing  wave  In  nzlayr 

dh2k1x=dh2xq1«k1xq2 
h2dk  1  xsh2xq1  -t-dkl  xq2 
h2dk2x-h2xq1'Klk2xq2 
dh2k2x=dh2xq1'»k2xq2 

call  xcadd(denax,cldqzl*cneg'*^dh2k1x,cldqzm«^cneg+h2dk1x) 
cal  I  xcadd(numax,cldqznW'h2dk2x,cldqzl'Hdh2k2x) 
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262 

263  c  calculate  contributions  to  cuaulants 

264 

265  8ui)x>cdlog(  'ql/dqi  jdz(  I  )‘^/dqi  jdz(  l«1 )  )'»2. 0d0*h2xq1 

266  call  xcadd(dhlx,dk2xq2,acoefx(ln1,M>^1xq2) 

267  dhlx«bcoefx(lMl,ia)'»dhlx 

268  int1x>2.0d0*dh2xq1-cldqzl 

269  call  xcadd(int1x,8unx,int1x) 

270  int2x>2.0d0*dhlx-cldqzni 

271  call  xcadd(suiix,int1x,int2x) 

272 

273  do  9030  l=nzlayr-1.2,-1 

274  lni1=l-1 

275  cldqzlscdlog(-dqi jdz(l)} 

276  cldqzRfxcdlogCdqi  jdz(lffll)) 

277 

278  c  calculate  q  and  associated  quantities  at  bottom  of  layer  I 

279 

280  qlscqi j(l,1)+zero*dqi j(l) 

281  call  Xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 

282  dk2xq1sdk2xq1't'cneg 

283  dk1xq1=dk1xq1-e13x 

284 

285  q2>cqi  j  ( Iml ,  2)-*’Zero*dqi  j  ( Iml ) 

286  call  Xcdai(-q2,k2xq2.dk2xq2.k1xq2,dk1xq2,h2xq2,dh2xq2) 

287  dk2xq2sdk2xq2'*’cneg 

288  dk1xq2>dk1xq2-e13x 

289  dh2xq2=dh2xq2-^e13x 

290 

291  c***** 

292  c  Calculate  acoefx(lm1,m),bcoefx(lm1,m)  and  cinulants 

293  c  using  continuity  relations  in  terms  of  the  linearly 

294  c  independent  functions  k1  and  k2 

295 

296  call  xcadd(hyx,k2xq1,acoefx(l,m)'«’k1xq1) 

297  call  xcadd(dhyx,dk2xq1,acoefx(l,ffl)4dk1xq1) 

298  k1dhyx>k1xq2^yx 

299  dk1hyx*dk1xq2+hyx 

300  dk2hyx=dk2xq2'*'hyx 

301  k2dhyxsk2xq2'»dhyx 

302 

303  call  xcadd(denax,cldqzl-cneg'»k1dhyx,cldqzm4’cneg'*'dk1hyx) 

304  call  xcadd(nunax,cldqznHdk2hyx,cldqzl'*'k2dhyx) 

305  acoefx(lm1,ffl)zn«jniax-denax 

306  call  xcadd(denbx,k2xq2,acoefx(lm1,m}+k1xq2) 
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307  nu«bx>bcoef  x(  I  ,«)'*'hyx 

308  dnuiibx>bcoefx(l,m}'Hjhyx 

309  bcoefx(  Ini  ,m)>nuRt>x-denbx 

310 

311  c  calculate  contributions  to  cuiulants. 

312 

313  int1x*cdlog(-q1/dqi  jdz(  I  )^/dqi  jclz(  Ini  ))4^2.0d0*nunbx 

314  call  xcadd(suiix,suBx,int1x) 

315  call  xcadd(dhlx,dk2xq2,acoefx(ln1,n)+dk1xq2) 

316  dhlx^bcoefxdnl.nl+dhlx 

317  int1xs2.0d0*dnijit>x-cldqzl 

318  int2xs2.0d0*dhlx-cldqzin 

319  call  xcadcKsunx.sunx.intlx) 

320  call  xcadd(stjnx,sijnx,int2x) 

321 

322  9030  continue 

323 

324  c***** 

325  c  if  t  equal  to  one  calculate  ground 

326  c  contribution  to  cuiulants  and  renormalize  bcoefx's 

327 

328  1=1 

329  q1=cqi j(l,1)>zero*dqi j(l) 

330  call  xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 

331  dk2xq1=dk2xq1«cneg 

332  dk1xq1=dk1xq1-e13x 

333 

334  call  xcadd(hyx,k2xq1,acoefx(l,in)'»k1xq1) 

335  call  xcadd(dhyx,dk2xq1,acoefx(l,ni)«dk1xq1) 

336  call  surf(q1,gannia,dganKlq} 

337  nunbx=bcoefx(l,ni)+hyx 

338  dnunbx=bcoefx(t,ni}'''dhyx 

339  int1x=cdlog(koaHav*d9affldq-q1/dqi  jdz(l))-'’2.0d0*nunt)x 

340  int2x=2.0d0*dnunbx-cdlog(-dqi  jdzd)) 

341  call  xcadd(sunx,sunx,int1x) 

342  call  xcadd(sunx,sumx,{nt2x} 

343 

344  c  renormalize  b's  so  that  height  gain  integrals  equal  unity. 

345 

346  rtsunx=.5d0*sunx 

347 

348  do  9020  11=1, nzlayr-1 

349  bcoefx(ll,ffl}=bcoefx(ll,m)*rtsunx 

350  9020  continue 

351 
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352  bcoefx(nzlayr,m}s-rtsuRU( 

353 

354 

355  return 

356  end 


« 


i 
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